!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine tvd_mmc (ro, mx, my, mz, bx, by, bz, en, az, po, gx, gy, x, y,                                                          &
   nx, ny, ng, n_neighs, dx, dy, dt, max_et, gm, cd, cfl, div_clean_flag, resistivity_flag,                                        &
   fxro_bnd, fxmx_bnd, fxmy_bnd, fxmz_bnd, fxbx_bnd, fxby_bnd, fxbz_bnd, fxen_bnd, fxpo_bnd,                                       &
   fyro_bnd, fymx_bnd, fymy_bnd, fymz_bnd, fybx_bnd, fyby_bnd, fybz_bnd, fyen_bnd, fypo_bnd)
!==================================================================================================================================|

   use interface_check_state
   implicit none

   integer(4), intent(in) :: nx, ny, ng, n_neighs, div_clean_flag, resistivity_flag
   real(8), intent(in) :: dx, dy, dt, gm, cd, cfl, x(nx), y(ny)
   real(8), intent(inout) :: max_et
   real(8), dimension(nx, ny), intent(inout) :: ro, mx, my, mz, bx, by, bz, en, az, po, gx, gy
   real(8), dimension(ny, n_neighs), intent(inout) ::                                                                              &
      fxro_bnd, fxmx_bnd, fxmy_bnd, fxmz_bnd, fxbx_bnd, fxby_bnd, fxbz_bnd, fxen_bnd, fxpo_bnd
   real(8), dimension(nx, n_neighs), intent(inout) ::                                                                              &
      fyro_bnd, fymx_bnd, fymy_bnd, fymz_bnd, fybx_bnd, fyby_bnd, fybz_bnd, fyen_bnd, fypo_bnd

   integer(4) :: i, j, ip1, im1, jp1, jm1

   real(8) :: vx, vy, vz, v2, b2, pr, gmm1, eps, c2, s2, ca2, cfx, cfy, max_speed
   real(8) :: ch, ch2, cp2
   real(8) :: dt_dx, dt_dy, dx2, dy2

   real(8), dimension(nx, ny) :: fxro, fxmx, fxmy, fxmz, fxbx, fxby, fxbz, fxen, fxpo
   real(8), dimension(nx, ny) :: fxro_p, fxmx_p, fxmy_p, fxmz_p, fxbx_p, fxby_p, fxbz_p, fxen_p, fxpo_p
   real(8), dimension(nx, ny) :: fxro_n, fxmx_n, fxmy_n, fxmz_n, fxbx_n, fxby_n, fxbz_n, fxen_n, fxpo_n
   real(8), dimension(nx, ny) :: fxro_up, fxmx_up, fxmy_up, fxmz_up, fxbx_up, fxby_up, fxbz_up, fxen_up, fxpo_up
   real(8), dimension(nx, ny) :: ro_pre, mx_pre, my_pre, mz_pre, bx_pre, by_pre, bz_pre, en_pre, az_pre, po_pre
   real(8), dimension(nx, ny) :: fyro, fymx, fymy, fymz, fybx, fyby, fybz, fyen, fypo
   real(8), dimension(nx, ny) :: fyro_p, fymx_p, fymy_p, fymz_p, fybx_p, fyby_p, fybz_p, fyen_p, fypo_p
   real(8), dimension(nx, ny) :: fyro_n, fymx_n, fymy_n, fymz_n, fybx_n, fyby_n, fybz_n, fyen_n, fypo_n
   real(8), dimension(nx, ny) :: fyro_up, fymx_up, fymy_up, fymz_up, fybx_up, fyby_up, fybz_up, fyen_up, fypo_up

   real(8), dimension(nx, ny) :: bx_resis, by_resis, bz_resis, en_resis, az_resis, sx, sy !, te_cond_x, te_cond_y

   real(8), dimension(nx, ny) :: fxro_pre_p, fxmx_pre_p, fxmy_pre_p, fxmz_pre_p,                                                   &
      fxbx_pre_p, fxby_pre_p, fxbz_pre_p, fxen_pre_p, fxpo_pre_p
   real(8), dimension(nx, ny) :: fxro_pre_n, fxmx_pre_n, fxmy_pre_n, fxmz_pre_n,                                                   &
      fxbx_pre_n, fxby_pre_n, fxbz_pre_n, fxen_pre_n, fxpo_pre_n
   real(8), dimension(nx, ny) :: fxro_lmt_p, fxmx_lmt_p, fxmy_lmt_p, fxmz_lmt_p,                                                   &
      fxbx_lmt_p, fxby_lmt_p, fxbz_lmt_p, fxen_lmt_p, fxpo_lmt_p
   real(8), dimension(nx, ny) :: fxro_lmt_n, fxmx_lmt_n, fxmy_lmt_n, fxmz_lmt_n,                                                   &
      fxbx_lmt_n, fxby_lmt_n, fxbz_lmt_n, fxen_lmt_n, fxpo_lmt_n
   real(8), dimension(nx, ny) :: fyro_pre_p, fymx_pre_p, fymy_pre_p, fymz_pre_p,                                                   &
      fybx_pre_p, fyby_pre_p, fybz_pre_p, fyen_pre_p, fypo_pre_p
   real(8), dimension(nx, ny) :: fyro_pre_n, fymx_pre_n, fymy_pre_n, fymz_pre_n,                                                   &
      fybx_pre_n, fyby_pre_n, fybz_pre_n, fyen_pre_n, fypo_pre_n
   real(8), dimension(nx, ny) :: fyro_lmt_p, fymx_lmt_p, fymy_lmt_p, fymz_lmt_p,                                                   &
      fybx_lmt_p, fyby_lmt_p, fybz_lmt_p, fyen_lmt_p, fypo_lmt_p
   real(8), dimension(nx, ny) :: fyro_lmt_n, fymx_lmt_n, fymy_lmt_n, fymz_lmt_n,                                                   &
      fybx_lmt_n, fyby_lmt_n, fybz_lmt_n, fyen_lmt_n, fypo_lmt_n

!----------------------------------------------------------------------------------------------------------------------------------|
   eps = 1d-12
   gmm1 = gm - 1.0d0
   dt_dx = dt / dx
   dt_dy = dt / dy
   dx2 = dx * 2.0d0
   dy2 = dy * 2.0d0
   ch = cfl * min (dx, dy) / dt * div_clean_flag
   ch2 = ch * ch
!   cp2 = ch * cd + eps
   cp2 = -dt * ch2 / log (cd) + eps

!   do j = 1, ny
!   do i = 1, nx
!      vx = mx(i, j) / ro(i, j) 
!      vy = my(i, j) / ro(i, j) 
!      vz = mz(i, j) / ro(i, j) 
!      b2 = bx(i, j) * bx(i, j) + by(i, j) * by(i, j) + bz(i, j) * bz(i, j)
!      v2 = vx * vx + vy * vy + vz * vz
!      pr = (en(i, j) - v2 * ro(i, j) / 2.0d0 - b2 / 2.0d0) * gmm1
!      c2 = gm * pr
!      s2 = c2 + b2
!      ca2 = bx(i, j) * bx(i, j)
!      cfx = sqrt((s2 + sqrt(s2 * s2 - 4.0d0 * c2 * ca2)) / ro(i, j)  / 2.0d0)
!      sx(i, j) = cfx + abs (vx)
!      ca2 = by(i, j) * by(i, j)
!      cfy = sqrt((s2 + sqrt(s2 * s2 - 4.0d0 * c2 * ca2)) / ro(i, j)  / 2.0d0)
!      sy(i, j) = cfy + abs (vy)
!      if (isnan(sx(i, j)) .or. sx(i, j) .gt. 1.0d0 / eps .or. sx(i, j) .lt. -1.0d0 / eps) sx(i, j) = 0.0d0
!      if (isnan(sy(i, j)) .or. sy(i, j) .gt. 1.0d0 / eps .or. sy(i, j) .lt. -1.0d0 / eps) sy(i, j) = 0.0d0
!   enddo
!   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  step one
!----------------------------------------------------------------------------------------------------------------------------------|
   do j = 1, ny
   do i = 1, nx
      vx = mx(i, j) / ro(i, j)
      vy = my(i, j) / ro(i, j)
      vz = mz(i, j) / ro(i, j)
      b2 = bx(i, j) * bx(i, j) + by(i, j) * by(i, j) + bz(i, j) * bz(i, j)
      v2 = vx * vx + vy * vy + vz * vz
      pr = (en(i, j) - v2 * ro(i, j) / 2.0d0 - b2 / 2.0d0) * gmm1
      c2 = gm * pr
      s2 = c2 + b2
      ca2 = bx(i, j) * bx(i, j)
      cfx = sqrt((s2 + sqrt(s2 * s2 - 4.0d0 * c2 * ca2)) / ro(i, j)  / 2.0d0)
      sx(i, j) = cfx + abs (vx)
      ca2 = by(i, j) * by(i, j)
      cfy = sqrt((s2 + sqrt(s2 * s2 - 4.0d0 * c2 * ca2)) / ro(i, j)  / 2.0d0)
      sy(i, j) = cfy + abs (vy)
      if (isnan(sx(i, j)) .or. sx(i, j) .gt. 1.0d0 / eps .or. sx(i, j) .lt. -1.0d0 / eps) sx(i, j) = 0.0d0
      if (isnan(sy(i, j)) .or. sy(i, j) .gt. 1.0d0 / eps .or. sy(i, j) .lt. -1.0d0 / eps) sy(i, j) = 0.0d0
   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!     flux
!----------------------------------------------------------------------------------------------------------------------------------|
   do j = 1, ny
   do i = 1, nx
      vx = mx(i, j) / ro(i, j)
      vy = my(i, j) / ro(i, j)
      vz = mz(i, j) / ro(i, j)
      b2 = bx(i, j) * bx(i, j) + by(i, j) * by(i, j) + bz(i, j) * bz(i, j)
      v2 = vx * vx + vy * vy + vz * vz
      pr = (en(i, j) - v2 * ro(i, j) / 2.0d0 - b2 / 2.0d0) * gmm1

      fxro(i, j) = mx(i, j)
      fxmx(i, j) = vx * mx(i, j) + pr + b2 / 2.0d0 - bx(i, j) * bx(i, j)
      fxmy(i, j) = vx * my(i, j) - bx(i, j) * by(i, j)
      fxmz(i, j) = vx * mz(i, j) - bx(i, j) * bz(i, j)
      fxbx(i, j) = po(i, j)
      fxby(i, j) = vx * by(i, j) - bx(i, j) * vy
      fxbz(i, j) = vx * bz(i, j) - bx(i, j) * vz
      fxen(i, j) = (en(i, j) + pr + b2 / 2.0d0) * vx - (bx(i, j) * vx + by(i, j) * vy + bz(i, j) * vz) * bx(i, j)
      fxpo(i, j) = ch2 * bx(i, j)

      fyro(i, j) = my(i, j)
      fymx(i, j) = vy * mx(i, j) - by(i, j) * bx(i, j)
      fymy(i, j) = vy * my(i, j) + pr + b2 / 2.0d0 - by(i, j) * by(i, j)
      fymz(i, j) = vy * mz(i, j) - by(i, j) * bz(i, j)
      fybx(i, j) = vy * bx(i, j) - by(i, j) * vx
      fyby(i, j) = po(i, j)
      fybz(i, j) = vy * bz(i, j) - by(i, j) * vz
      fyen(i, j) = (en(i, j) + pr + b2 / 2.0d0) * vy - (bx(i, j) * vx + by(i, j) * vy + bz(i, j) * vz) * by(i, j)
      fypo(i, j) = ch2 * by(i, j)

!----------------------------------------------------------------------------------------------------------------------------------|
!     flux split
!----------------------------------------------------------------------------------------------------------------------------------|
      if (i .eq. 1) then
         max_speed = max (sx(i, j), sx(i + 1, j))
      elseif (i .eq. nx) then
         max_speed = max (sx(i - 1, j), sx(i, j))
      else
         max_speed = maxval (sx(i - 1 : i + 1, j))
      endif

      fxro_p(i, j) = (fxro(i, j) + max_speed * ro(i, j)) * 0.5d0
      fxro_n(i, j) = (fxro(i, j) - max_speed * ro(i, j)) * 0.5d0
      fxmx_p(i, j) = (fxmx(i, j) + max_speed * mx(i, j)) * 0.5d0
      fxmx_n(i, j) = (fxmx(i, j) - max_speed * mx(i, j)) * 0.5d0
      fxmy_p(i, j) = (fxmy(i, j) + max_speed * my(i, j)) * 0.5d0
      fxmy_n(i, j) = (fxmy(i, j) - max_speed * my(i, j)) * 0.5d0
      fxmz_p(i, j) = (fxmz(i, j) + max_speed * mz(i, j)) * 0.5d0
      fxmz_n(i, j) = (fxmz(i, j) - max_speed * mz(i, j)) * 0.5d0
      fxbx_p(i, j) = (fxbx(i, j) + max_speed * bx(i, j)) * 0.5d0
      fxbx_n(i, j) = (fxbx(i, j) - max_speed * bx(i, j)) * 0.5d0
      fxby_p(i, j) = (fxby(i, j) + max_speed * by(i, j)) * 0.5d0
      fxby_n(i, j) = (fxby(i, j) - max_speed * by(i, j)) * 0.5d0
      fxbz_p(i, j) = (fxbz(i, j) + max_speed * bz(i, j)) * 0.5d0
      fxbz_n(i, j) = (fxbz(i, j) - max_speed * bz(i, j)) * 0.5d0
      fxen_p(i, j) = (fxen(i, j) + max_speed * en(i, j)) * 0.5d0
      fxen_n(i, j) = (fxen(i, j) - max_speed * en(i, j)) * 0.5d0
      fxpo_p(i, j) = (fxpo(i, j) + max_speed * po(i, j)) * 0.5d0
      fxpo_n(i, j) = (fxpo(i, j) - max_speed * po(i, j)) * 0.5d0

      if (j .eq. 1) then
         max_speed = max (sy(i, j), sy(i, j + 1))
      elseif (j .eq. ny) then
         max_speed = max (sy(i, j - 1), sy(i, j))
      else
         max_speed = maxval (sy(i, j - 1 : j + 1))
      endif

      fyro_p(i, j) = (fyro(i, j) + max_speed * ro(i, j)) * 0.5d0
      fyro_n(i, j) = (fyro(i, j) - max_speed * ro(i, j)) * 0.5d0
      fymx_p(i, j) = (fymx(i, j) + max_speed * mx(i, j)) * 0.5d0
      fymx_n(i, j) = (fymx(i, j) - max_speed * mx(i, j)) * 0.5d0
      fymy_p(i, j) = (fymy(i, j) + max_speed * my(i, j)) * 0.5d0
      fymy_n(i, j) = (fymy(i, j) - max_speed * my(i, j)) * 0.5d0
      fymz_p(i, j) = (fymz(i, j) + max_speed * mz(i, j)) * 0.5d0
      fymz_n(i, j) = (fymz(i, j) - max_speed * mz(i, j)) * 0.5d0
      fybx_p(i, j) = (fybx(i, j) + max_speed * bx(i, j)) * 0.5d0
      fybx_n(i, j) = (fybx(i, j) - max_speed * bx(i, j)) * 0.5d0
      fyby_p(i, j) = (fyby(i, j) + max_speed * by(i, j)) * 0.5d0
      fyby_n(i, j) = (fyby(i, j) - max_speed * by(i, j)) * 0.5d0
      fybz_p(i, j) = (fybz(i, j) + max_speed * bz(i, j)) * 0.5d0
      fybz_n(i, j) = (fybz(i, j) - max_speed * bz(i, j)) * 0.5d0
      fyen_p(i, j) = (fyen(i, j) + max_speed * en(i, j)) * 0.5d0
      fyen_n(i, j) = (fyen(i, j) - max_speed * en(i, j)) * 0.5d0
      fypo_p(i, j) = (fypo(i, j) + max_speed * po(i, j)) * 0.5d0
      fypo_n(i, j) = (fypo(i, j) - max_speed * po(i, j)) * 0.5d0
   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  first-order upwind numerical flux
!----------------------------------------------------------------------------------------------------------------------------------|
   do j = 1, ny - 1
      jp1 = j + 1
      do i = 1, nx - 1
         ip1 = i + 1
         fxro_up(i, j) = fxro_p(i, j) + fxro_n(ip1, j)
         fxmx_up(i, j) = fxmx_p(i, j) + fxmx_n(ip1, j)
         fxmy_up(i, j) = fxmy_p(i, j) + fxmy_n(ip1, j)
         fxmz_up(i, j) = fxmz_p(i, j) + fxmz_n(ip1, j)
         fxbx_up(i, j) = fxbx_p(i, j) + fxbx_n(ip1, j)
         fxby_up(i, j) = fxby_p(i, j) + fxby_n(ip1, j)
         fxbz_up(i, j) = fxbz_p(i, j) + fxbz_n(ip1, j)
         fxen_up(i, j) = fxen_p(i, j) + fxen_n(ip1, j)
         fxpo_up(i, j) = fxpo_p(i, j) + fxpo_n(ip1, j)

         fyro_up(i, j) = fyro_p(i, j) + fyro_n(i, jp1)
         fymx_up(i, j) = fymx_p(i, j) + fymx_n(i, jp1)
         fymy_up(i, j) = fymy_p(i, j) + fymy_n(i, jp1)
         fymz_up(i, j) = fymz_p(i, j) + fymz_n(i, jp1)
         fybx_up(i, j) = fybx_p(i, j) + fybx_n(i, jp1)
         fyby_up(i, j) = fyby_p(i, j) + fyby_n(i, jp1)
         fybz_up(i, j) = fybz_p(i, j) + fybz_n(i, jp1)
         fyen_up(i, j) = fyen_p(i, j) + fyen_n(i, jp1)
         fypo_up(i, j) = fypo_p(i, j) + fypo_n(i, jp1)
      enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  predictor-corrector step
!----------------------------------------------------------------------------------------------------------------------------------|
   do j = 2, ny
      jm1 = j - 1
      do i = 2, nx
         im1 = i - 1
         ro_pre(i, j) = ro(i, j) - dt_dx * (fxro_up(i, j) - fxro_up(im1, j)) - dt_dy * (fyro_up(i, j) - fyro_up(i, jm1))
         mx_pre(i, j) = mx(i, j) - dt_dx * (fxmx_up(i, j) - fxmx_up(im1, j)) - dt_dy * (fymx_up(i, j) - fymx_up(i, jm1))
         my_pre(i, j) = my(i, j) - dt_dx * (fxmy_up(i, j) - fxmy_up(im1, j)) - dt_dy * (fymy_up(i, j) - fymy_up(i, jm1))
         mz_pre(i, j) = mz(i, j) - dt_dx * (fxmz_up(i, j) - fxmz_up(im1, j)) - dt_dy * (fymz_up(i, j) - fymz_up(i, jm1))
         bx_pre(i, j) = bx(i, j) - dt_dx * (fxbx_up(i, j) - fxbx_up(im1, j)) - dt_dy * (fybx_up(i, j) - fybx_up(i, jm1))
         by_pre(i, j) = by(i, j) - dt_dx * (fxby_up(i, j) - fxby_up(im1, j)) - dt_dy * (fyby_up(i, j) - fyby_up(i, jm1))
         bz_pre(i, j) = bz(i, j) - dt_dx * (fxbz_up(i, j) - fxbz_up(im1, j)) - dt_dy * (fybz_up(i, j) - fybz_up(i, jm1))
         en_pre(i, j) = en(i, j) - dt_dx * (fxen_up(i, j) - fxen_up(im1, j)) - dt_dy * (fyen_up(i, j) - fyen_up(i, jm1))
         po_pre(i, j) = po(i, j) - dt_dx * (fxpo_up(i, j) - fxpo_up(im1, j)) - dt_dy * (fypo_up(i, j) - fypo_up(i, jm1))
         az_pre(i, j) = az(i, j)
      enddo
   enddo

   call check_state (ro_pre, mx_pre, my_pre, mz_pre, bx_pre, by_pre, bz_pre, en_pre, x, y, nx, ny, gm)

   if (resistivity_flag .eq. 0) then
      do j = 1, ny
      do i = 1, nx
         bx_resis(i, j) = 0.0d0
         by_resis(i, j) = 0.0d0
         bz_resis(i, j) = 0.0d0
         en_resis(i, j) = 0.0d0
         az_resis(i, j) = 0.0d0
      enddo
      enddo
      max_et = 0.0d0
   else
      call resistivity (bx_resis, by_resis, bz_resis, en_resis, az_resis,                                                          &
         ro, mx, my, mz, bx, by, bz, en, x, y, dx, dy, nx, ny, max_et)
   endif

!   call thermal_conduction (te_cond_x, te_cond_y, ro, mx, my, mz, bx, by, bz, en, thermal_amp, gm, dx, dy, nx, ny)

   call source_term (mx_pre, my_pre, mz_pre, bx_pre, by_pre, bz_pre, en_pre, az_pre, po_pre,                                       &
      ro, mx, my, bx, by, bz, po, gx, gy, bx_resis, by_resis, bz_resis, en_resis, az_resis,                                        &
      ch2, cp2, dx, dy, dt, nx, ny, div_clean_flag)

!----------------------------------------------------------------------------------------------------------------------------------|
!  step two
!----------------------------------------------------------------------------------------------------------------------------------|
   do j = 1, ny
   do i = 1, nx
      vx = mx_pre(i, j) / ro_pre(i, j)
      vy = my_pre(i, j) / ro_pre(i, j)
      vz = mz_pre(i, j) / ro_pre(i, j)
      b2 = bx_pre(i, j) * bx_pre(i, j) + by_pre(i, j) * by_pre(i, j) + bz_pre(i, j) * bz_pre(i, j)
      v2 = vx * vx + vy * vy + vz * vz
      pr = (en_pre(i, j) - v2 * ro_pre(i, j) / 2.0d0 - b2 / 2.0d0) * gmm1
      c2 = gm * pr
      s2 = c2 + b2
      ca2 = bx_pre(i, j) * bx_pre(i, j)
      cfx = sqrt((s2 + sqrt(s2 * s2 - 4.0d0 * c2 * ca2)) / ro_pre(i, j)  / 2.0d0)
      sx(i, j) = cfx + abs (vx)
      ca2 = by(i, j) * by(i, j)
      cfy = sqrt((s2 + sqrt(s2 * s2 - 4.0d0 * c2 * ca2)) / ro_pre(i, j)  / 2.0d0)
      sy(i, j) = cfy + abs (vy)
      if (isnan(sx(i, j)) .or. sx(i, j) .gt. 1.0d0 / eps .or. sx(i, j) .lt. -1.0d0 / eps) sx(i, j) = 0.0d0
      if (isnan(sy(i, j)) .or. sy(i, j) .gt. 1.0d0 / eps .or. sy(i, j) .lt. -1.0d0 / eps) sy(i, j) = 0.0d0
   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!     flux
!----------------------------------------------------------------------------------------------------------------------------------|
   do j = 1, ny
   do i = 1, nx
      vx = mx_pre(i, j) / ro_pre(i, j)
      vy = my_pre(i, j) / ro_pre(i, j)
      vz = mz_pre(i, j) / ro_pre(i, j)
      b2 = bx_pre(i, j) * bx_pre(i, j) + by_pre(i, j) * by_pre(i, j) + bz_pre(i, j) * bz_pre(i, j)
      v2 = vx * vx + vy * vy + vz * vz
      pr = (en_pre(i, j) - v2 * ro_pre(i, j) / 2.0d0 - b2 / 2.0d0) * gmm1

      fxro(i, j) = mx_pre(i, j)
      fxmx(i, j) = vx * mx_pre(i, j) + pr + b2 / 2.0d0 - bx_pre(i, j) * bx_pre(i, j)
      fxmy(i, j) = vx * my_pre(i, j) - bx_pre(i, j) * by_pre(i, j)
      fxmz(i, j) = vx * mz_pre(i, j) - bx_pre(i, j) * bz_pre(i, j)
      fxbx(i, j) = po_pre(i, j)
      fxby(i, j) = vx * by_pre(i, j) - bx_pre(i, j) * vy
      fxbz(i, j) = vx * bz_pre(i, j) - bx_pre(i, j) * vz
      fxen(i, j) = (en_pre(i, j) + pr + b2 / 2) * vx - (bx_pre(i, j) * vx + by_pre(i, j) * vy + bz_pre(i, j) * vz) * bx_pre(i, j)
      fxpo(i, j) = ch2 * bx_pre(i, j)

      fyro(i, j) = my_pre(i, j)
      fymx(i, j) = vy * mx_pre(i, j) - by_pre(i, j) * bx_pre(i, j)
      fymy(i, j) = vy * my_pre(i, j) + pr + b2 / 2.0d0 - by_pre(i, j) * by_pre(i, j)
      fymz(i, j) = vy * mz_pre(i, j) - by_pre(i, j) * bz_pre(i, j)
      fybx(i, j) = vy * bx_pre(i, j) - by_pre(i, j) * vx
      fyby(i, j) = po_pre(i, j)
      fybz(i, j) = vy * bz_pre(i, j) - by_pre(i, j) * vz
      fyen(i, j) = (en_pre(i, j) + pr + b2 / 2) * vy - (bx_pre(i, j) * vx + by_pre(i, j) * vy + bz_pre(i, j) * vz) * by_pre(i, j)
      fypo(i, j) = ch2 * by_pre(i, j)

!----------------------------------------------------------------------------------------------------------------------------------|
!     flux split
!----------------------------------------------------------------------------------------------------------------------------------|
      if (i .eq. 1) then
         max_speed = max (sx(i, j), sx(i + 1, j))
      elseif (i .eq. nx) then
         max_speed = max (sx(i - 1, j), sx(i, j))
      else
         max_speed = maxval (sx(i - 1 : i + 1, j))
      endif

      fxro_pre_p(i, j) = (fxro(i, j) + max_speed * ro_pre(i, j)) * 0.5d0
      fxro_pre_n(i, j) = (fxro(i, j) - max_speed * ro_pre(i, j)) * 0.5d0
      fxmx_pre_p(i, j) = (fxmx(i, j) + max_speed * mx_pre(i, j)) * 0.5d0
      fxmx_pre_n(i, j) = (fxmx(i, j) - max_speed * mx_pre(i, j)) * 0.5d0
      fxmy_pre_p(i, j) = (fxmy(i, j) + max_speed * my_pre(i, j)) * 0.5d0
      fxmy_pre_n(i, j) = (fxmy(i, j) - max_speed * my_pre(i, j)) * 0.5d0
      fxmz_pre_p(i, j) = (fxmz(i, j) + max_speed * mz_pre(i, j)) * 0.5d0
      fxmz_pre_n(i, j) = (fxmz(i, j) - max_speed * mz_pre(i, j)) * 0.5d0
      fxbx_pre_p(i, j) = (fxbx(i, j) + max_speed * bx_pre(i, j)) * 0.5d0
      fxbx_pre_n(i, j) = (fxbx(i, j) - max_speed * bx_pre(i, j)) * 0.5d0
      fxby_pre_p(i, j) = (fxby(i, j) + max_speed * by_pre(i, j)) * 0.5d0
      fxby_pre_n(i, j) = (fxby(i, j) - max_speed * by_pre(i, j)) * 0.5d0
      fxbz_pre_p(i, j) = (fxbz(i, j) + max_speed * bz_pre(i, j)) * 0.5d0
      fxbz_pre_n(i, j) = (fxbz(i, j) - max_speed * bz_pre(i, j)) * 0.5d0
      fxen_pre_p(i, j) = (fxen(i, j) + max_speed * en_pre(i, j)) * 0.5d0
      fxen_pre_n(i, j) = (fxen(i, j) - max_speed * en_pre(i, j)) * 0.5d0
      fxpo_pre_p(i, j) = (fxpo(i, j) + max_speed * po_pre(i, j)) * 0.5d0
      fxpo_pre_n(i, j) = (fxpo(i, j) - max_speed * po_pre(i, j)) * 0.5d0

      if (j .eq. 1) then
         max_speed = max (sy(i, j), sy(i, j + 1))
      elseif (j .eq. ny) then
         max_speed = max (sy(i, j - 1), sy(i, j))
      else
         max_speed = maxval (sy(i, j - 1 : j + 1))
      endif
      fyro_pre_p(i, j) = (fyro(i, j) + max_speed * ro_pre(i, j)) * 0.5d0
      fyro_pre_n(i, j) = (fyro(i, j) - max_speed * ro_pre(i, j)) * 0.5d0
      fymx_pre_p(i, j) = (fymx(i, j) + max_speed * mx_pre(i, j)) * 0.5d0
      fymx_pre_n(i, j) = (fymx(i, j) - max_speed * mx_pre(i, j)) * 0.5d0
      fymy_pre_p(i, j) = (fymy(i, j) + max_speed * my_pre(i, j)) * 0.5d0
      fymy_pre_n(i, j) = (fymy(i, j) - max_speed * my_pre(i, j)) * 0.5d0
      fymz_pre_p(i, j) = (fymz(i, j) + max_speed * mz_pre(i, j)) * 0.5d0
      fymz_pre_n(i, j) = (fymz(i, j) - max_speed * mz_pre(i, j)) * 0.5d0
      fybx_pre_p(i, j) = (fybx(i, j) + max_speed * bx_pre(i, j)) * 0.5d0
      fybx_pre_n(i, j) = (fybx(i, j) - max_speed * bx_pre(i, j)) * 0.5d0
      fyby_pre_p(i, j) = (fyby(i, j) + max_speed * by_pre(i, j)) * 0.5d0
      fyby_pre_n(i, j) = (fyby(i, j) - max_speed * by_pre(i, j)) * 0.5d0
      fybz_pre_p(i, j) = (fybz(i, j) + max_speed * bz_pre(i, j)) * 0.5d0
      fybz_pre_n(i, j) = (fybz(i, j) - max_speed * bz_pre(i, j)) * 0.5d0
      fyen_pre_p(i, j) = (fyen(i, j) + max_speed * en_pre(i, j)) * 0.5d0
      fyen_pre_n(i, j) = (fyen(i, j) - max_speed * en_pre(i, j)) * 0.5d0
      fypo_pre_p(i, j) = (fypo(i, j) + max_speed * po_pre(i, j)) * 0.5d0
      fypo_pre_n(i, j) = (fypo(i, j) - max_speed * po_pre(i, j)) * 0.5d0
   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  antidiffusive flux term
!----------------------------------------------------------------------------------------------------------------------------------|
   do j = 1, ny - 1
      jp1 = j + 1
      do i = 1, nx - 1
         ip1 = i + 1
         fxro_p(i, j) = fxro_pre_p(ip1, j) - fxro_p(i, j)
         fxro_n(i, j) = fxro_n(ip1, j) - fxro_pre_n(i, j)
         fxmx_p(i, j) = fxmx_pre_p(ip1, j) - fxmx_p(i, j)
         fxmx_n(i, j) = fxmx_n(ip1, j) - fxmx_pre_n(i, j)
         fxmy_p(i, j) = fxmy_pre_p(ip1, j) - fxmy_p(i, j)
         fxmy_n(i, j) = fxmy_n(ip1, j) - fxmy_pre_n(i, j)
         fxmz_p(i, j) = fxmz_pre_p(ip1, j) - fxmz_p(i, j)
         fxmz_n(i, j) = fxmz_n(ip1, j) - fxmz_pre_n(i, j)
         fxbx_p(i, j) = fxbx_pre_p(ip1, j) - fxbx_p(i, j)
         fxbx_n(i, j) = fxbx_n(ip1, j) - fxbx_pre_n(i, j)
         fxby_p(i, j) = fxby_pre_p(ip1, j) - fxby_p(i, j)
         fxby_n(i, j) = fxby_n(ip1, j) - fxby_pre_n(i, j)
         fxbz_p(i, j) = fxbz_pre_p(ip1, j) - fxbz_p(i, j)
         fxbz_n(i, j) = fxbz_n(ip1, j) - fxbz_pre_n(i, j)
         fxen_p(i, j) = fxen_pre_p(ip1, j) - fxen_p(i, j)
         fxen_n(i, j) = fxen_n(ip1, j) - fxen_pre_n(i, j)
         fxpo_p(i, j) = fxpo_pre_p(ip1, j) - fxpo_p(i, j)
         fxpo_n(i, j) = fxpo_n(ip1, j) - fxpo_pre_n(i, j)

         fyro_p(i, j) = fyro_pre_p(i, jp1) - fyro_p(i, j)
         fyro_n(i, j) = fyro_n(i, jp1) - fyro_pre_n(i, j)
         fymx_p(i, j) = fymx_pre_p(i, jp1) - fymx_p(i, j)
         fymx_n(i, j) = fymx_n(i, jp1) - fymx_pre_n(i, j)
         fymy_p(i, j) = fymy_pre_p(i, jp1) - fymy_p(i, j)
         fymy_n(i, j) = fymy_n(i, jp1) - fymy_pre_n(i, j)
         fymz_p(i, j) = fymz_pre_p(i, jp1) - fymz_p(i, j)
         fymz_n(i, j) = fymz_n(i, jp1) - fymz_pre_n(i, j)
         fybx_p(i, j) = fybx_pre_p(i, jp1) - fybx_p(i, j)
         fybx_n(i, j) = fybx_n(i, jp1) - fybx_pre_n(i, j)
         fyby_p(i, j) = fyby_pre_p(i, jp1) - fyby_p(i, j)
         fyby_n(i, j) = fyby_n(i, jp1) - fyby_pre_n(i, j)
         fybz_p(i, j) = fybz_pre_p(i, jp1) - fybz_p(i, j)
         fybz_n(i, j) = fybz_n(i, jp1) - fybz_pre_n(i, j)
         fyen_p(i, j) = fyen_pre_p(i, jp1) - fyen_p(i, j)
         fyen_n(i, j) = fyen_n(i, jp1) - fyen_pre_n(i, j)
         fypo_p(i, j) = fypo_pre_p(i, jp1) - fypo_p(i, j)
         fypo_n(i, j) = fypo_n(i, jp1) - fypo_pre_n(i, j)
      enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  non-linear flux-limiter
!----------------------------------------------------------------------------------------------------------------------------------|
   do j = 2, ny
      jm1 = j - 1
      do i = 2, nx
         im1 = i - 1
         if (abs(fxro_p(i, j)) .eq. 0.0d0) fxro_p(i, j) = eps
         if (abs(fxmx_p(i, j)) .eq. 0.0d0) fxmx_p(i, j) = eps
         if (abs(fxmy_p(i, j)) .eq. 0.0d0) fxmy_p(i, j) = eps
         if (abs(fxmz_p(i, j)) .eq. 0.0d0) fxmz_p(i, j) = eps
         if (abs(fxbx_p(i, j)) .eq. 0.0d0) fxbx_p(i, j) = eps
         if (abs(fxby_p(i, j)) .eq. 0.0d0) fxby_p(i, j) = eps
         if (abs(fxbz_p(i, j)) .eq. 0.0d0) fxbz_p(i, j) = eps
         if (abs(fxen_p(i, j)) .eq. 0.0d0) fxen_p(i, j) = eps
         if (abs(fxpo_p(i, j)) .eq. 0.0d0) fxpo_p(i, j) = eps

         fxro_lmt_p(i, j) = fxro_p(im1, j) / fxro_p(i, j)
         fxmx_lmt_p(i, j) = fxmx_p(im1, j) / fxmx_p(i, j)
         fxmy_lmt_p(i, j) = fxmy_p(im1, j) / fxmy_p(i, j)
         fxmz_lmt_p(i, j) = fxmz_p(im1, j) / fxmz_p(i, j)
         fxbx_lmt_p(i, j) = fxbx_p(im1, j) / fxbx_p(i, j)
         fxby_lmt_p(i, j) = fxby_p(im1, j) / fxby_p(i, j)
         fxbz_lmt_p(i, j) = fxbz_p(im1, j) / fxbz_p(i, j)
         fxen_lmt_p(i, j) = fxen_p(im1, j) / fxen_p(i, j)
         fxpo_lmt_p(i, j) = fxpo_p(im1, j) / fxpo_p(i, j)

         fxro_lmt_p(i, j) = min(1.d0, abs(fxro_lmt_p(i, j)))
         fxmx_lmt_p(i, j) = min(1.d0, abs(fxmx_lmt_p(i, j)))
         fxmy_lmt_p(i, j) = min(1.d0, abs(fxmy_lmt_p(i, j)))
         fxmz_lmt_p(i, j) = min(1.d0, abs(fxmz_lmt_p(i, j)))
         fxbx_lmt_p(i, j) = min(1.d0, abs(fxbx_lmt_p(i, j)))
         fxby_lmt_p(i, j) = min(1.d0, abs(fxby_lmt_p(i, j)))
         fxbz_lmt_p(i, j) = min(1.d0, abs(fxbz_lmt_p(i, j)))
         fxen_lmt_p(i, j) = min(1.d0, abs(fxen_lmt_p(i, j)))
         fxpo_lmt_p(i, j) = min(1.d0, abs(fxpo_lmt_p(i, j)))

         if (abs(fyro_p(i, j)) .eq. 0.0d0) fyro_p(i, j) = eps
         if (abs(fymx_p(i, j)) .eq. 0.0d0) fymx_p(i, j) = eps
         if (abs(fymy_p(i, j)) .eq. 0.0d0) fymy_p(i, j) = eps
         if (abs(fymz_p(i, j)) .eq. 0.0d0) fymz_p(i, j) = eps
         if (abs(fybx_p(i, j)) .eq. 0.0d0) fybx_p(i, j) = eps
         if (abs(fyby_p(i, j)) .eq. 0.0d0) fyby_p(i, j) = eps
         if (abs(fybz_p(i, j)) .eq. 0.0d0) fybz_p(i, j) = eps
         if (abs(fyen_p(i, j)) .eq. 0.0d0) fyen_p(i, j) = eps
         if (abs(fypo_p(i, j)) .eq. 0.0d0) fypo_p(i, j) = eps

         fyro_lmt_p(i, j) = fyro_p(i, jm1) / fyro_p(i, j)
         fymx_lmt_p(i, j) = fymx_p(i, jm1) / fymx_p(i, j)
         fymy_lmt_p(i, j) = fymy_p(i, jm1) / fymy_p(i, j)
         fymz_lmt_p(i, j) = fymz_p(i, jm1) / fymz_p(i, j)
         fybx_lmt_p(i, j) = fybx_p(i, jm1) / fybx_p(i, j)
         fyby_lmt_p(i, j) = fyby_p(i, jm1) / fyby_p(i, j)
         fybz_lmt_p(i, j) = fybz_p(i, jm1) / fybz_p(i, j)
         fyen_lmt_p(i, j) = fyen_p(i, jm1) / fyen_p(i, j)
         fypo_lmt_p(i, j) = fypo_p(i, jm1) / fypo_p(i, j)

         fyro_lmt_p(i, j) = min(1.d0, abs(fyro_lmt_p(i, j)))
         fymx_lmt_p(i, j) = min(1.d0, abs(fymx_lmt_p(i, j)))
         fymy_lmt_p(i, j) = min(1.d0, abs(fymy_lmt_p(i, j)))
         fymz_lmt_p(i, j) = min(1.d0, abs(fymz_lmt_p(i, j)))
         fybx_lmt_p(i, j) = min(1.d0, abs(fybx_lmt_p(i, j)))
         fyby_lmt_p(i, j) = min(1.d0, abs(fyby_lmt_p(i, j)))
         fybz_lmt_p(i, j) = min(1.d0, abs(fybz_lmt_p(i, j)))
         fyen_lmt_p(i, j) = min(1.d0, abs(fyen_lmt_p(i, j)))
         fypo_lmt_p(i, j) = min(1.d0, abs(fypo_lmt_p(i, j)))
      enddo
   enddo

   do j = 1, ny - 1
      jp1 = j + 1
      do i = 1, nx - 1
         ip1 = i + 1
         if (abs(fxro_n(i, j)) .eq. 0.0d0) fxro_n(i, j) = eps
         if (abs(fxmx_n(i, j)) .eq. 0.0d0) fxmx_n(i, j) = eps
         if (abs(fxmy_n(i, j)) .eq. 0.0d0) fxmy_n(i, j) = eps
         if (abs(fxmz_n(i, j)) .eq. 0.0d0) fxmz_n(i, j) = eps
         if (abs(fxbx_n(i, j)) .eq. 0.0d0) fxbx_n(i, j) = eps
         if (abs(fxby_n(i, j)) .eq. 0.0d0) fxby_n(i, j) = eps
         if (abs(fxbz_n(i, j)) .eq. 0.0d0) fxbz_n(i, j) = eps
         if (abs(fxen_n(i, j)) .eq. 0.0d0) fxen_n(i, j) = eps
         if (abs(fxpo_n(i, j)) .eq. 0.0d0) fxpo_n(i, j) = eps

         fxro_lmt_n(i, j) = fxro_n(ip1, j) / fxro_n(i, j)
         fxmx_lmt_n(i, j) = fxmx_n(ip1, j) / fxmx_n(i, j)
         fxmy_lmt_n(i, j) = fxmy_n(ip1, j) / fxmy_n(i, j)
         fxmz_lmt_n(i, j) = fxmz_n(ip1, j) / fxmz_n(i, j)
         fxbx_lmt_n(i, j) = fxbx_n(ip1, j) / fxbx_n(i, j)
         fxby_lmt_n(i, j) = fxby_n(ip1, j) / fxby_n(i, j)
         fxbz_lmt_n(i, j) = fxbz_n(ip1, j) / fxbz_n(i, j)
         fxen_lmt_n(i, j) = fxen_n(ip1, j) / fxen_n(i, j)
         fxpo_lmt_n(i, j) = fxpo_n(ip1, j) / fxpo_n(i, j)

         fxro_lmt_n(i, j) = min(1.d0, abs(fxro_lmt_n(i, j)))
         fxmx_lmt_n(i, j) = min(1.d0, abs(fxmx_lmt_n(i, j)))
         fxmy_lmt_n(i, j) = min(1.d0, abs(fxmy_lmt_n(i, j)))
         fxmz_lmt_n(i, j) = min(1.d0, abs(fxmz_lmt_n(i, j)))
         fxbx_lmt_n(i, j) = min(1.d0, abs(fxbx_lmt_n(i, j)))
         fxby_lmt_n(i, j) = min(1.d0, abs(fxby_lmt_n(i, j)))
         fxbz_lmt_n(i, j) = min(1.d0, abs(fxbz_lmt_n(i, j)))
         fxen_lmt_n(i, j) = min(1.d0, abs(fxen_lmt_n(i, j)))
         fxpo_lmt_n(i, j) = min(1.d0, abs(fxpo_lmt_n(i, j)))

         if (abs(fyro_n(i, j)) .eq. 0.0d0) fyro_n(i, j) = eps
         if (abs(fymx_n(i, j)) .eq. 0.0d0) fymx_n(i, j) = eps
         if (abs(fymy_n(i, j)) .eq. 0.0d0) fymy_n(i, j) = eps
         if (abs(fymz_n(i, j)) .eq. 0.0d0) fymz_n(i, j) = eps
         if (abs(fybx_n(i, j)) .eq. 0.0d0) fybx_n(i, j) = eps
         if (abs(fyby_n(i, j)) .eq. 0.0d0) fyby_n(i, j) = eps
         if (abs(fybz_n(i, j)) .eq. 0.0d0) fybz_n(i, j) = eps
         if (abs(fyen_n(i, j)) .eq. 0.0d0) fyen_n(i, j) = eps
         if (abs(fypo_n(i, j)) .eq. 0.0d0) fypo_n(i, j) = eps

         fyro_lmt_n(i, j) = fyro_n(i, jp1) / fyro_n(i, j)
         fymx_lmt_n(i, j) = fymx_n(i, jp1) / fymx_n(i, j)
         fymy_lmt_n(i, j) = fymy_n(i, jp1) / fymy_n(i, j)
         fymz_lmt_n(i, j) = fymz_n(i, jp1) / fymz_n(i, j)
         fybx_lmt_n(i, j) = fybx_n(i, jp1) / fybx_n(i, j)
         fyby_lmt_n(i, j) = fyby_n(i, jp1) / fyby_n(i, j)
         fybz_lmt_n(i, j) = fybz_n(i, jp1) / fybz_n(i, j)
         fyen_lmt_n(i, j) = fyen_n(i, jp1) / fyen_n(i, j)
         fypo_lmt_n(i, j) = fypo_n(i, jp1) / fypo_n(i, j)

         fyro_lmt_n(i, j) = min(1.d0, abs(fyro_lmt_n(i, j)))
         fymx_lmt_n(i, j) = min(1.d0, abs(fymx_lmt_n(i, j)))
         fymy_lmt_n(i, j) = min(1.d0, abs(fymy_lmt_n(i, j)))
         fymz_lmt_n(i, j) = min(1.d0, abs(fymz_lmt_n(i, j)))
         fybx_lmt_n(i, j) = min(1.d0, abs(fybx_lmt_n(i, j)))
         fyby_lmt_n(i, j) = min(1.d0, abs(fyby_lmt_n(i, j)))
         fybz_lmt_n(i, j) = min(1.d0, abs(fybz_lmt_n(i, j)))
         fyen_lmt_n(i, j) = min(1.d0, abs(fyen_lmt_n(i, j)))
         fypo_lmt_n(i, j) = min(1.d0, abs(fypo_lmt_n(i, j)))
      enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!   second-order numerical flux
!----------------------------------------------------------------------------------------------------------------------------------|
   do j = 1, ny
   do i = 1, nx
      fxro(i, j) = fxro_up(i, j) + (fxro_lmt_p(i, j) * fxro_p(i, j) - fxro_lmt_n(i, j) * fxro_n(i, j)) * 0.5d0
      fxmx(i, j) = fxmx_up(i, j) + (fxmx_lmt_p(i, j) * fxmx_p(i, j) - fxmx_lmt_n(i, j) * fxmx_n(i, j)) * 0.5d0
      fxmy(i, j) = fxmy_up(i, j) + (fxmy_lmt_p(i, j) * fxmy_p(i, j) - fxmy_lmt_n(i, j) * fxmy_n(i, j)) * 0.5d0
      fxmz(i, j) = fxmz_up(i, j) + (fxmz_lmt_p(i, j) * fxmz_p(i, j) - fxmz_lmt_n(i, j) * fxmz_n(i, j)) * 0.5d0
      fxbx(i, j) = fxbx_up(i, j) + (fxbx_lmt_p(i, j) * fxbx_p(i, j) - fxbx_lmt_n(i, j) * fxbx_n(i, j)) * 0.5d0
      fxby(i, j) = fxby_up(i, j) + (fxby_lmt_p(i, j) * fxby_p(i, j) - fxby_lmt_n(i, j) * fxby_n(i, j)) * 0.5d0
      fxbz(i, j) = fxbz_up(i, j) + (fxbz_lmt_p(i, j) * fxbz_p(i, j) - fxbz_lmt_n(i, j) * fxbz_n(i, j)) * 0.5d0
      fxen(i, j) = fxen_up(i, j) + (fxen_lmt_p(i, j) * fxen_p(i, j) - fxen_lmt_n(i, j) * fxen_n(i, j)) * 0.5d0
      fxpo(i, j) = fxpo_up(i, j) + (fxpo_lmt_p(i, j) * fxpo_p(i, j) - fxpo_lmt_n(i, j) * fxpo_n(i, j)) * 0.5d0

      fyro(i, j) = fyro_up(i, j) + (fyro_lmt_p(i, j) * fyro_p(i, j) - fyro_lmt_n(i, j) * fyro_n(i, j)) * 0.5d0
      fymx(i, j) = fymx_up(i, j) + (fymx_lmt_p(i, j) * fymx_p(i, j) - fymx_lmt_n(i, j) * fymx_n(i, j)) * 0.5d0
      fymy(i, j) = fymy_up(i, j) + (fymy_lmt_p(i, j) * fymy_p(i, j) - fymy_lmt_n(i, j) * fymy_n(i, j)) * 0.5d0
      fymz(i, j) = fymz_up(i, j) + (fymz_lmt_p(i, j) * fymz_p(i, j) - fymz_lmt_n(i, j) * fymz_n(i, j)) * 0.5d0
      fybx(i, j) = fybx_up(i, j) + (fybx_lmt_p(i, j) * fybx_p(i, j) - fybx_lmt_n(i, j) * fybx_n(i, j)) * 0.5d0
      fyby(i, j) = fyby_up(i, j) + (fyby_lmt_p(i, j) * fyby_p(i, j) - fyby_lmt_n(i, j) * fyby_n(i, j)) * 0.5d0
      fybz(i, j) = fybz_up(i, j) + (fybz_lmt_p(i, j) * fybz_p(i, j) - fybz_lmt_n(i, j) * fybz_n(i, j)) * 0.5d0
      fyen(i, j) = fyen_up(i, j) + (fyen_lmt_p(i, j) * fyen_p(i, j) - fyen_lmt_n(i, j) * fyen_n(i, j)) * 0.5d0
      fypo(i, j) = fypo_up(i, j) + (fypo_lmt_p(i, j) * fypo_p(i, j) - fypo_lmt_n(i, j) * fypo_n(i, j)) * 0.5d0
   enddo
   enddo

   do j = 2, ny
      jm1 = j - 1
      do i = 2, nx
         im1 = i - 1
         ro(i, j) = ro(i, j) - dt_dx * (fxro(i, j) - fxro(im1, j)) - dt_dy * (fyro(i, j) - fyro(i, jm1))
         mx(i, j) = mx(i, j) - dt_dx * (fxmx(i, j) - fxmx(im1, j)) - dt_dy * (fymx(i, j) - fymx(i, jm1))
         my(i, j) = my(i, j) - dt_dx * (fxmy(i, j) - fxmy(im1, j)) - dt_dy * (fymy(i, j) - fymy(i, jm1))
         mz(i, j) = mz(i, j) - dt_dx * (fxmz(i, j) - fxmz(im1, j)) - dt_dy * (fymz(i, j) - fymz(i, jm1))
         bx(i, j) = bx(i, j) - dt_dx * (fxbx(i, j) - fxbx(im1, j)) - dt_dy * (fybx(i, j) - fybx(i, jm1))
         by(i, j) = by(i, j) - dt_dx * (fxby(i, j) - fxby(im1, j)) - dt_dy * (fyby(i, j) - fyby(i, jm1))
         bz(i, j) = bz(i, j) - dt_dx * (fxbz(i, j) - fxbz(im1, j)) - dt_dy * (fybz(i, j) - fybz(i, jm1))
         en(i, j) = en(i, j) - dt_dx * (fxen(i, j) - fxen(im1, j)) - dt_dy * (fyen(i, j) - fyen(i, jm1))
         po(i, j) = po(i, j) - dt_dx * (fxpo(i, j) - fxpo(im1, j)) - dt_dy * (fypo(i, j) - fypo(i, jm1))
      enddo
   enddo

   if (resistivity_flag .eq. 0) then
      do j = 1, ny
      do i = 1, nx
         bx_resis(i, j) = 0.0d0
         by_resis(i, j) = 0.0d0
         bz_resis(i, j) = 0.0d0
         en_resis(i, j) = 0.0d0
         az_resis(i, j) = 0.0d0
      enddo
      enddo
      max_et = 0.0d0
   else
      call resistivity (bx_resis, by_resis, bz_resis, en_resis, az_resis,                                                          &
         ro_pre, mx_pre, my_pre, mz_pre, bx_pre, by_pre, bz_pre, en_pre, x, y, dx, dy, nx, ny, max_et)
   endif

!   call thermal_conduction (te_cond_x, te_cond_y, ro_pre, mx_pre, my_pre, mz_pre, bx_pre, by_pre, bz_pre, en_pre,                  &
!      thermal_amp, gm, dx, dy, nx, ny)

   call source_term (mx, my, mz, bx, by, bz, en, az, po,                                                                           &
      ro_pre, mx_pre, my_pre, bx_pre, by_pre, bz_pre, po_pre, gx, gy,                                                              &
      bx_resis, by_resis, bz_resis, en_resis, az_resis,                                                                            &
      ch2, cp2, dx, dy, dt, nx, ny, div_clean_flag)

!----------------------------------------------------------------------------------------------------------------------------------|
!  bnd value
!----------------------------------------------------------------------------------------------------------------------------------|
   do j = 1, ny
      fxro_bnd(j, 1) = fxro(ng, j)
      fxro_bnd(j, 2) = fxro(nx - ng, j)

      fxmx_bnd(j, 1) = fxmx(ng, j)
      fxmx_bnd(j, 2) = fxmx(nx - ng, j)

      fxmy_bnd(j, 1) = fxmy(ng, j)
      fxmy_bnd(j, 2) = fxmy(nx - ng, j)

      fxmz_bnd(j, 1) = fxmz(ng, j)
      fxmz_bnd(j, 2) = fxmz(nx - ng, j)

      fxbx_bnd(j, 1) = fxbx(ng, j)
      fxbx_bnd(j, 2) = fxbx(nx - ng, j)

      fxby_bnd(j, 1) = fxby(ng, j)
      fxby_bnd(j, 2) = fxby(nx - ng, j)

      fxbz_bnd(j, 1) = fxbz(ng, j)
      fxbz_bnd(j, 2) = fxbz(nx - ng, j)

      fxen_bnd(j, 1) = fxen(ng, j)
      fxen_bnd(j, 2) = fxen(nx - ng, j)

      fxpo_bnd(j, 1) = fxpo(ng, j)
      fxpo_bnd(j, 2) = fxpo(nx - ng, j)
   enddo

   do i = 1, nx
      fyro_bnd(i, 1) = fyro(i, ng)
      fyro_bnd(i, 2) = fyro(i, ny - ng)

      fymx_bnd(i, 1) = fymx(i, ng)
      fymx_bnd(i, 2) = fymx(i, ny - ng)

      fymy_bnd(i, 1) = fymy(i, ng)
      fymy_bnd(i, 2) = fymy(i, ny - ng)

      fymz_bnd(i, 1) = fymz(i, ng)
      fymz_bnd(i, 2) = fymz(i, ny - ng)

      fybx_bnd(i, 1) = fybx(i, ng)
      fybx_bnd(i, 2) = fybx(i, ny - ng)

      fyby_bnd(i, 1) = fyby(i, ng)
      fyby_bnd(i, 2) = fyby(i, ny - ng)

      fybz_bnd(i, 1) = fybz(i, ng)
      fybz_bnd(i, 2) = fybz(i, ny - ng)

      fyen_bnd(i, 1) = fyen(i, ng)
      fyen_bnd(i, 2) = fyen(i, ny - ng)

      fypo_bnd(i, 1) = fypo(i, ng)
      fypo_bnd(i, 2) = fypo(i, ny - ng)
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine tvd_mmc
